Journal of Mathematical Biology

This paper concerns the inverse problem of characterising the state of a bioreactor from observations. In laboratory settings, the bioreactor is represented by a device called a chemostat. We consider a differential description of the evolution of the state of the chemostat under environmental fluctuations. First, we model the state evolution as a stochastic process driven by Brownian motion. Under this model, our best knowledge about the state of the chemostat is described by its probability distribution in time, given the distribution of the initial state. The corresponding probability density function solves a deterministic partial differential equation (PDE), the Kolmogorov forward equation. While this provides a probabilistic description, incorporating an observation process allows for a more refined characterisation of the state. More formally, we are interested in obtaining the distribution of the state conditional on an observation process as the solution to a filtering problem, with the corresponding conditional probability density function solving a non-linear stochastic PDE, the Kushner-Stratonovich equation. This paper focuses on the pathwise formulation of this filtering problem in which inferences about the state are obtained conditional on a fixed stream of observations. We establish the existence and uniqueness of solutions to the governing differential equations, ensuring well-posedness before presenting numerical approximations. We approximate the pathwise solution to the filtering problem by combining the finite difference and splitting methods for solving PDEs, and then compare the approximated solution with results from a linearisation method and a classical sequential Monte Carlo method.

Since its introduction by Briat, Gupta and Khammash, the antithetic feedback controller design has attracted considerable attention in both theoretical and experimental systems biology. The case in which the plant is a two-dimensional linear system (making the closed-loop system a nonlinear four-dimensional system) has been analyzed in much detail. This system has a unique equilibrium but, depending on parameters, it may exhibit periodic orbits. An interesting open question is whether other dynamical behaviors, such as chaotic attractors, might be possible for some parameter choices. This note shows that, for any parameter choices, every bounded trajectory satisfies a Poincare-Bendixson property. The analysis is based on the recently introduced notion of k-cooperative dynamical systems. It is shown that the model is a strongly 2-cooperative system, implying that the dynamics in the omega-limit set of any precompact solution is conjugate to the dynamics in a compact invariant subset of a two-dimensional Lipschitz dynamical system, thus precluding chaotic and other strange attractors.

A fresh approach to the dynamics of gene assemblies is presented. Central to the exposition are the concepts of: high value genes; correlated activity; and the orderly unfolding of gene dynamics; and especially dynamic mode decomposition, DMD, a remarkable new tool for dissecting dynamics. This program is carried out, in detail, for the Orlando et al yeast database (Orlando et al. 2008). It is shown that the yeast cell division cycle, CDC, requires no more than a six dimensional space, formed by three complex temporal modal pairs, each associated with characteristic aspects of the cell cycle: (1) A mother cell cohort that follows a fast clock; (2) A daughter cell cohort that follows a slower clock; (3) inherent gene expression, unrelated to the CDC. A derived set of sixty high-value genes serves as a model for the correlated unfolding of gene activity. Confirmation of our results comes from an independent database, and other considerations. The present analysis, leads naturally, to a Fourier description, for the sparsely sampled data. From this, resolved peak times of gene expression are obtained. This in turn leads to prediction of precise times of expression in the unfolding of the CDC genes. The activation of each gene appears as uncoupled dynamics from the mother and daughter cohorts, of different durations. These deliberations lead to detailed estimates of the fraction of mother and daughter cells, specific estimates of their maturation periods, and specific estimates of the number of genes in these cells. An algorithmic framework for yeast modeling is proposed, and based on the new analyses, a range of theoretical ideas and new experiments are suggested. A Supplement contains additional material and other perspectives.

The exact mechanisms behind many morphogenic processes are still a mystery. Mechanical cues, such as curvature, play an important role when tissue or cell shape is formed. In this work, we derive and analyze a mechanochemical model. This particular spatially one-dimensional model describes the deformation of a tissue- or cell surface over time, which is driven by a morphogen that locally induces curvature. The model consists of two PDEs with periodic boundary conditions; one reaction-diffusion equation for the morphogen and one PDE that describes the dynamics of the curve, derived by taking the L2-gradient flow of the Helfrich energy. We analyze the possible steady states of this model using geometric singular perturbation theory. It turns out that the strength of interaction between the morphogen and the curvature plays a key role in the type of possible steady state solutions. In the case of weak interaction, the geometry of the slow manifolds allows only for (in space) slowly changing periodic orbits that lay completely on one slow manifold. In the case of strong interaction, there exist multiple front solutions: periodic orbits that jump between different slow manifolds. The singular skeletons of the steady state solutions do not meet the required consistency conditions for the curvature, a priori indicating that the solutions might not be observable. The observability and stability are investigated further using numerical simulation.

It is well known that the presence of an incoherent feedforward loop (IFFL) in a network may give rise to a steady state non-monotonic dose response. This note shows that the converse implication does not hold. It gives an example of a three-dimensional system that has no IFFLs, yet its dose response is bell-shaped. It also studies under what conditions the result is true for two-dimensional systems, in the process recovering, in far more generality, a result given in the T-cell activation literature.

The analysis of large gene and metabolic networks is often hindered by unknown biochemical parameters and the nonlinear nature of classical S-system models. To address this, we introduce a framework based on combinatorial toric geometry computed with tools such as Normaliz, SageMath, it is worth mentioning this technique in not restrictive to integer vectors, there exists a natural extension to real geometries. Unlike traditional approaches, which rely on parameter dependent fixed points, our method constructs a Topological Environment derived from the dual space of kinetic orders, leading to what we call orthogonal enzyme kinetics. Within this topological setting, fixed points are computed on the algebraic torus, enabling the transformation of nonlinear dynamics into linear forms. Importantly, these fixed points are independent of kinetic parameters and depend only on network topology and interaction signs. Applying this methodology to gene circuits involved in circadian rhythms, we reproduce previously reported oscillatory physiologies.

This note describes the outcome of an epidemic in a heterogeneous population with a very simple structure. The population is split into two sections in which the epidemic runs a different course. The reproduction numbers of the two epidemics are unobserved, only the overall reproduction number is known. For such a population the outcome of the epidemic can be as well far worse as far better than expected on the basis of the overall reproduction number. By considering a very simple model of this population some calculations are feasible under general assumptions on the epidemic itself. These calculations show in which direction models based on the overall reproduction number can misrepresent the real-world situation.

We discuss some critical events of the origins of life using a mathematical model and simulation studies. We find that for a replicating population of RNA molecules participating in template-directed polymerization, the hitting and establishment of a high-fidelity replicator depends critically on the polymerase fitness and sequence specificity landscapes and on genome dimension. Probability of hitting is dominated by polymerase landscape curvature, whereas hitting time is dominated by genome dimension. Surface chemistries, compartmentalization, and decay increase hitting times. These results suggest replication to be the first privileged function marking the start of Darwinian evolution, possibly in conjunction with clay minerals or preceded by metabolism, whose dynamics evolved mostly during the final period of the search.

The tetrapod limb skeleton is initiated in unpatterned limb bud mesenchyme by the formation of precartilage condensations. Here, based on time-lapse videographic analysis of a forming condensation in a high-density culture of chicken limb bud mesenchyme, we observe a phase transition to a more fluidized state for cells within spatial compacted foci (protocondensations that will progress to condensations), as reflected in their spatial confinement, cell-substratum interaction and speed of motion. Previous work showed that galectin-8 and galectin-1A, two proteins of the galactoside-binding galectin family, are the earliest determinants of this process in the chicken limb bud, and that their interactions in forming skeletogenic patterns of condensations can be interpreted mathematically through a reaction-diffusion-adhesion framework. Based on this framework, we use an ordinary differential equation-based approach to analyze the core switching modality of the galectin reaction network and characterize the states of the network independent of the diffusive and adhesive arms of the patterning mechanism. We identify two steady states where the concentrations of both galectins are respectively, negligible, and very high. An explicit Lyapunov function shows that there are no periodic solutions. For sigmoidal galectin production terms, the model exhibits a bistable switch that arises from a monostable state via saddle-node bifurcation. Our model therefore predicts that the galectin network exists in low and high expression states separated in space or time without any intermediate states. This provides a causal basis for the observed outside vs. inside transition observed in the in vitro video data. We performed a quantitative analysis of the distribution of galectin-1A in cultures of condensing chick limb mesenchymal cells and found that the interior of the protocondensations had concentrations of this protein (compared to the immediate exterior) over and above that expected from its higher cell density, consistent with the models predictions. The galectin-based patterning network is thus suggested, on theoretical grounds, to incorporate a core switch independent of any spatial or temporal dynamics, that drives the chondrogenic cell state transition in limb skeletogenesis.

This paper discusses the mean duration of a closed epidemic modeled by a discrete-time Markov chain. We develop a methodology for the efficient computation of the quantity of interest. The Markov chain model in consideration is bivariate, and is formally handled. We derive explicit terms for the probability to transition from one state to another, and prove that the chain is absorbing. The computation of the mean duration is translated to the computation of the expected hitting times to the set of absorbing states. We use the theory of absorbing Markov chains to derive a matrix formulation that gives way to an efficient algorithm to solve for the expected hitting times. This approach is instantiated in the form of a concrete algorithm, which is further optimized by using dynamic programming. Finally, we have implemented the method and tested it against the use of simulations to estimate mean durations.

A survey of the literature reveals notable discrepancies among the purported exact results for the spectra of stochastic gene expression models. For self-repressing gene circuits, previous studies ([Phys. Rev. Lett. 99, 108103 (2007)], [Phys. Rev. E 83,062902 (2011)], [J. Chem. Phys. 160, 074105 (2024)], and [bioRxiv 2025.02.05.635946 (2025)]) have provided different exact solutions for the eigenvalues of the generator matrix. In this work, we propose a unified Hilbert space framework for the spectral theory of stochastic gene expression. Based on this framework, we analytically derive the spectra for models of constitutive, bursty, and autoregulated gene expression. The eigenvalues and eigenvectors obtained are then used to construct an exact spectral representation of the time-dependent distribution of gene product numbers. The spectral gap between the zero eigenvalue and the first nonzero eigenvalue, which reflects the relaxation rate of the system towards its steady state, is then compared with the prediction of the deterministic model, and we find that deterministic modeling fails to capture the relaxation rate when autoregulation is strong. In particular, our results demonstrate that for infinite-dimensional operators such as in stochastic gene expression models, many conclusions in linear algebra do not apply, and one must rely on the modern theory of functional analysis.

Surprisingly, the discrete-time version of the general 1927 Kermack-McKendrick epidemic model has, to our knowledge, not been formulated in the literature, and we rectify this omission here. The discrete time version is as general and flexible as its continuous-time counterpart, and contains numerous compartmental models as special cases. In contrast to the continuous time version, the discrete time version of the model is very easy to implement computationally, and thus promises to become a powerful tool for exploring control scenarios for specific infectious diseases. To demonstrate the potential, we investigate numerically how the incidence-peak size depends on model ingredients. We find that, with the same reproduction number and initial speed of epidemic spread, compartmental models systematically predict lower peak sizes than models that use a fixed duration for the latent and infectious periods.

Melanoma is a cancer of the melanocyte, known to have an ability to readily switch between different transcriptional cell states that convey different phenotypic properties (e.g. hyper-differentiated, neural crest-like). This ability is believed to underpin intratumour heterogeneity and plastic adaptation, which contributes to resistance to therapy and immune evasion of the tumour. Therefore, understanding the mechanisms underlying acquisition of transcriptional cell states and cell-state switching is crucial for the development of therapies. We model a minimal gene regulatory network comprising three key transcription factors, whose varying gene expression encodes different melanoma cell states, and use deterministic spatiotemporal differential-equation models to study gene-expression dynamics. We exploit an approximation, based on cooperative binding of transcription factors, in which the models are piecewise-linear. We classify stable states of the local model in a biologically relevant manner and, using a naive model of intercellular communication, we explore how a population of cells can take on a shared characteristic through travelling waves of gene expression. We derive a condition determining which characteristic will become dominant, under sufficiently strong cell-cell signalling, which creates a partition of parameter space.

AO_SCPLOWBSTRACTC_SCPLOWCompartmental models have long served as important tools in mathematical epidemiology, with their usefulness highlighted by the recent COVID-19 pandemic. However, most of the classical models fail to account for certain features of this disease and others like it, such as the ability of exposed individuals to recover without becoming infectious, or the possibility that asymptomatic individuals can indeed transmit the disease but at a lesser rate than the symptomatic. Furthermore, the rise of new disease variants and the imperfection of vaccines suggest that concept of endemic equilibrium is perhaps more pertinent than that of herd immunity. Here we propose a new compartmental epidemiological model and study its equilibria, characterizing the stability of both the endemic and disease-free equilibria in terms of the basic reproductive number. Moreover, we introduce a second compartmental model, generalizing our first, which accounts for vaccinated individuals, and begin an analysis of its equilibria.

In a variety of practical applications, there is a need to investigate diffusion or reaction-diffusion processes on complex structures, including brain networks, that can be modeled as weighted undirected and directed graphs. As an instance, the celebrated Fisher-Kolmogorov-Petrovsky-Piskunov (Fisher-KPP) reaction-diffusion equation are becoming increasingly popular for use in graph frameworks by substituting the standard graph Laplacian operator for the continuous one to study the progression of neurodegenerative diseases such as tauopathies including Alzheimers disease (AD). However, due to the porous structure of neuronal fibers, the spreading of toxic species can be governed by an anomalous diffusion process rather than a normal one, and if this is the case, the standard graph Laplacian cannot adequately describe the dynamics of the spreading process. To capture such more complicated dynamics, we propose a diffusion equation with a nonlinear Laplacian operator and a generalization of the Fisher-KPP reaction-diffusion equation on undirected and directed networks using extensions of fractional polynomial (FP) functions. A complete analysis is also provided for the extended FP diffusion equation, including existence, uniqueness, and convergence of solutions, as well as stability of equilibria. Moreover, for the extended FP Fisher-KPP reaction-diffusion equation, we derive a family of positively invariant sets allowing us to establish existence, uniqueness, and boundedness of solutions. Finally, we conclude by investigating nonlinear diffusion on a directed one-dimensional lattice and then modeling tauopathy progression in the mouse brain to gain a deeper understanding of the potential applications of the proposed extended FP equations.

Matrix games under time constraints are natural extensions of matrix games. They consider the fact that, in addition to the payoff, a pairwise interaction has a further consequence for the contestants. Namely, both players have to wait for some time before becoming fit to participate in a subsequent interaction. Every matrix game can be assigned a continuous dynamical system (the replicator equation) which describes how the frequencies of different phenotypes evolve in the population. One of the fundamental theorems of evolutionary matrix games asserts that the state corresponding to an evolutionarily stable strategy is an asymptotically stable rest point of the replicator equation (Taylor and Yonker 1978, Hofbauer et al. 1979, Zeeman 1980). Garay et al. (2018) and Varga et al. (2020) generalized the statement to two-strategy and, in some particular cases, three- or more strategy matrix games under time constraints. However, the question of whether the implication holds in general remained open. Here examples are provided demonstrating that the answer is no. Moreover, we point out through the rock-scissor-paper game that arbitrary small differences between waiting times can destabilize the rest point corresponding to an ESS. It is also shown that a stable limit cycle can arise around the unstable rest point in a supercritical Hopf bifurcation. Mathematics Subject Classification91A22, 92D15, 92D25, 91A80, 91A05, 91A10, 91A40, 92D40

Gene regulatory networks (GRN) control the expression levels of proteins in cells, and understanding their dynamics is key to potentially controlling disease processes. Steady states of GRNs are interpreted as cellular phenotypes, and the first step in understanding GRN dynamics is describing the collection of steady states the network can support in different conditions. We consider a collection of all monotone Boolean function models compatible with a given GRN, and ask which steady states are supported by most models. We find that for networks with no negative loops, there is an explicit hierarchy in the prevalence of individual steady states, as well as the prevalence of bistability and multistability. The key insight that we use is that monotone Boolean models supporting a given equilibrium are a product of prime ideals and prime filters of the lattices of monotone Boolean functions. To illustrate our result, we show that in the EMT network associated with cancer metastasis, the most common equilibria correspond to epithelial (E) and mesenchymal (M) states, and the bistability between them is the most common bistability among all network-compatible monotone Boolean models. Author summaryCells adjust their behavior in response to external inputs via networks of genes that regulate each others expression, until they arrive at a new steady state. Each interacting network of genes can behave in different ways that depend on internal and external cellular conditions. In this paper we consider, for a given network, an entire collection of particular type of models (monotone Boolean models) that represent all different ways that network can behave. Then, for any given state a network can potentially be in, we describe all monotone Boolean models that have that state as a steady state. We consider those states that are supported by more models to more likely represent the states that the network will be in. We apply our approach to EMT network that is important in cancer metastasis. We show that the most common steady states are those correspond to epithelial and mesenchymal states, and that the bistability between these two states is the most common bistability. This confirms the experimental results that these are the most common states of the EMT network.

We study the mechanism and effects of ideal free distributions (IFDs) on an ecological community consisting of n prey and m predator species, for any positive integers n and m, by considering the corresponding diffusive Lotka-Volterra system with time-periodic coefficients. We define a notion of joint IFD in a timeperiodic environment, and give necessary and sufficient conditions for it to be achieved by a subcollection of prey and predator species with suitable dispersal strategies. Next, we show, via construction of a Lyapunov function, that such dispersal strategies are evolutionarily stable, in the sense that if a subcollection of prey and predator species adopts an ideal free dispersal strategy, then the total community must converge to an IFD for large time; if a unique combination of prey-predator species adopts an ideal free strategy, then it can drive all other species to extinction. Conversely, if a combination of prey-predator species adopts a non-ideal free dispersal strategy, then it can be invaded by some suitable mutant strategies. Our results provide insight into the evolution of spatial distribution of ecological communities with predator-prey interactions.

Oncolytic virotherapy is a promising targeted cancer treatment that employs viruses, which selectively infect tumor cells. Although its clinical efficiency has remained limited and it is often used in conjunction with other therapies, advances in genetic engineering have produced stronger and more selective viral strains, prompting continued interest in their dynamics. In particular, previous studies have noted that viruses with sufficiently high replication rates can induce oscillations reminiscent of predator-prey systems. Here, we extend this analysis to the spatial domain by starting from an established tumor-virus reaction-diffusion model, performing a center-manifold reduction that incorporates nonlinear terms to derive a complex Ginzburg-Landau amplitude equation, and estimating its parameters directly from the original kinetics. This reduced normal form equation explains the emergence of experimentally observed patterns -- such as hollow rings and target waves -- and shows that, at longer timescales, these patterns naturally evolve toward spiral waves and a turbulent regime. Our work provides a mechanistic link between the kinetic Hopf bifurcation and the rich spatiotemporal structures observed in oncolytic virotherapy models, suggesting that these patterns are not numerical artifacts but an intrinsic feature of the system.

AO_SCPLOWBSTRACTC_SCPLOWWe investigate which evolutionary histories can potentially be reconstructed from sufficiently long DNA sequences by studying the identifiability of phylogenetic networks from sequence data generated under site independent models of molecular evolution. While previous work in the field has established the identifiability of phylogenetic trees and level-1 networks, networks with non-overlapping reticulation cycles, less is known about more complex network structures. In this work, we extend identifiability results to network classes that include pairs of tangled reticulations. Our main result shows that binary semi-directed level-2 phylogenetic networks are generically identifiable under the Jukes-Cantor model, provided they are triangle-free and strongly tree-child. We also strengthen existing identifiability results for level-1 networks, showing that the number of reticulation nodes is generically identifiable under the Jukes-Cantor model. In addition, we present more general identifiability results that do not restrict the network level at all and hold for the Jukes-Cantor as well as for the Kimura-2-Parameter model. Specifically, we demonstrate that any two binary semi-directed networks that display different sets of 4-leaf subtrees (quartets) are distinguishable. This has direct implications for the identifiability of a networks reticulated components (blobs). We show that the tree-of-blobs of a network, the global branching structure of the network, is identifiable, as well as the circular ordering of the subnetworks around each blob, for networks in which edges do not cross and taxa are on the outside.